data= load('iris_i.dat' );
data=data';
[dim N ]= size(data);
c= input( 'enter c (# of Classes) in Fuzzy_c_means algo :');
m= input( 'enter m ,the exponential weight :');
epsilon= input( 'enter epsilon, the termination criterion :');
flag= 1;
Normalize= [];
U0= rand(c,N);
Normalize= sum(U0,1);
for j= 1:N,
U0(:,j)= U0(:,j)/ Normalize(j);
end % end for

U1= zeros(c,N);
v= [];
for i= 1 :c,
    v(i).mean= zeros(1,dim);
    v(i).denom= 1;
end % end for
% computation of new v(i)'s
while flag
 U0= U0.^m;
 for i= 1:c,
     S= sum(U0,2);
     v(i).denom= S(i);
     numerator= zeros(1,dim);
     for j= 1 : N,
         numerator= numerator + (U0(i,j)*data(:,j))';
     end % end for j
     if v(i).denom ~= 0,
         v(i).mean= numerator/(v(i).denom);
     end % end if
 end % end for i
 
 % computation for U1
 inv_norm= zeros(c,N);
 inv_pow= 1/(m-1);
 for i= 1:c,
     for j= 1: N,
         temp= (sum(data(:,j)- v(i).mean').^2);
         if temp ~= 0,
             temp= temp^inv_pow;
         end % end if
         if temp ~= 0,
             inv_norm(i,j)= 1/temp;
         end %end if
     end % end for j
 end % end for i
 
 %while flag 
 for j=1:N,
     denom= 0;
     for k= 1:c,
         denom= denom + inv_norm(k,j);
     end % end for k
     for i= 1:c, 
         if denom ~= 0,    
             U1(i,j)= inv_norm(i,j)/denom ;
         end % end if
     end % end for i
 end %end for j
 
 if max(max( abs(U1- U0)))<= epsilon 
     flag= 0;
 end % end if
U0= U1;
end % end while

D=[];
C=[];
[Y,I]= max(U0,[],1);
for j= 1:N,
    D(j)=j;
    C(j)=I(j);
end
plot(D,C,'o');
